function [ P, Error ] = q4( )
%Q4 Summary of this function goes here
%   Implementation of MC method using Standard, B-B, and PCA constructions

% Parameters
S0 = 100; K = 100; r = 0.05; sigma = 0.3;
M = 10000; N = 32;
T = 1;
dt = 1/N;
disc = exp(-r*T);

% Standard Construction
Wt = randn(M,N);
factor = exp((r-0.5*sigma^2)*dt + sigma*sqrt(dt)*Wt);
St1 = S0*cumprod(factor,2);

P1 = disc*mean(max(mean(St1,2)-K,0));
Var1 = var(max(mean(St1,2)-K,0));

% Brownian Bridge Construction
Tn = 0:dt:1;
Zn = randn(M,N+1);
Wn = zeros(M,N+1);
m = log2(N);
h = 2^m; j_max = 1;
Wn(:,h+1) = sqrt(T)*Zn(:,h+1);
for k=1:m
    i_min = h/2; i = i_min + 1;
    l = 1; u = h+1;
    for j = 1:j_max
        a = ((Tn(u) - Tn(i))*Wn(:,l) + (Tn(i) - Tn(l))*Wn(:,u))/...
            (Tn(u) - Tn(l));
        b = sqrt((Tn(i) - Tn(l))*(Tn(u) - Tn(i))/(Tn(u) - Tn(l)));
        Wn(:,i) = a + b*Zn(:,i);
        i = i + h; l = l + h; u = u + h;
    end
    j_max = 2*j_max;
    h = i_min;
end
    
% Now Price based on Brownian Bridge
[x,y] = meshgrid(Tn,ones(M,1));
factor = (r - 0.5*sigma^2)*x + sigma*Wn; 
St2 = S0 * exp(factor(:,2:end));

P2 = disc*mean(max(mean(St2,2)-K,0));
Var2 = var(max(mean(St2,2)-K,0));


% Now the PCA construction

% Generate Covariance Matrix for Wt
[x,y] = meshgrid(Tn(2:end));
S = tril(x) + tril(x,-1)';
[V,D] = eig(S);
A = V*sqrt(D);
Z = randn(N,M);
W = A*Z;

% Now Price based on PCA path
[x,y] = meshgrid(Tn(2:end),ones(M,1));
factor = (r - 0.5*sigma^2)*x + sigma*W';
St3 = S0*exp(factor);

P3 = disc*mean(max(mean(St3,2)-K,0));
Var3 = var(max(mean(St3,2)-K,0));

P = [P1,P2,P3];
Error = sqrt([Var1, Var2, Var3]);

end

